#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "mpi.h"
//mpicc -std=c99 -o cmake-build/mpi/jacobi mpi/jacobi.c && mpirun -np 4 ./cmake-build/mpi/jacobi

#define matrixN (10)

double matrix[matrixN][matrixN];
double jacobi[matrixN][matrixN];

int main(int argc,char*argv[])
{
    int rank,size;
    MPI_Init(&argc,&argv);
    MPI_Comm_size(MPI_COMM_WORLD,&size);
    MPI_Comm_rank(MPI_COMM_WORLD,&rank);
    
    //进程0初始化输入矩阵和结果矩阵
    if(rank==0)
    {
        for(int i=0;i<matrixN;++i)
            for(int j=0;j<matrixN;j+=2)
                matrix[i][j]=i+j+1;
        
        for(int i=0;i<matrixN;++i)
            for(int j=0;j<matrixN;++j)
                jacobi[i][j]=0;
        
        for(int i=0;i<matrixN;++i,putchar('\n'))
            for(int j=0;j<matrixN;++j)
                printf("%6.2lf ",matrix[i][j]);
    }
    //向所有进程广播
    MPI_Bcast(matrix[0],matrixN*matrixN,MPI_DOUBLE,0,MPI_COMM_WORLD);
    MPI_Bcast(jacobi[0],matrixN*matrixN,MPI_DOUBLE,0,MPI_COMM_WORLD);
    
    //迭代
    for(int i=rank;i<matrixN;i+=size)
    {
        printf("Thread %d process row %d\n",rank,i);
        if(i==0||i==matrixN-1)
        {
            for(int j=0;j<matrixN;++j)
                jacobi[i][j]=matrix[i][j];
        }
        else
        {
            jacobi[i][0]=matrix[i][0];
            for(int j=1;j<matrixN-1;++j)
                jacobi[i][j]=(matrix[i-1][j]+matrix[i+1][j]+matrix[i][j-1]+matrix[i][j+1])/4;
            jacobi[i][matrixN-1]=matrix[i][matrixN-1];
        }
        fflush(stdout);
    }
    //发回结果
    MPI_Barrier(MPI_COMM_WORLD);
    MPI_Allreduce(MPI_IN_PLACE,jacobi[0],matrixN*matrixN,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
    
    //输出结果
    fflush(stdout);
    if(rank==0)
    {
        for(int i=0;i<matrixN;++i,putchar('\n'))
            for(int j=0;j<matrixN;++j)
                printf("%6.2lf ",jacobi[i][j]);
    }
    
    MPI_Finalize();
    return 0;
}
